Influence of solvent quality on effective pair potentials between polymers in solution 
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Solutions of interacting linear polymers are mapped onto a system of "soft" spherical particles 
interacting via an effective pair potential. This coarse-graining reduces the individual monomer-level 
description to a problem involving only the centers of mass (CM) of the polymer coils. The effective 
pair potentials are derived by inverting the CM pair distribution function, generated in Monte Carlo 
simulations, using the hypernetted chain (HNC) closure. The method, previously devised for the self- 
avoiding walk model of polymers in good solvent, is extended to the case of polymers in solvents of 
variable quality by adding a finite nearest-neighbor monomer-monomer attraction to the previous 
model and varying the temperature. The resulting effective pair potential is found to depend 
strongly on temperature and polymer concentration. At low concentration the effective interaction 
becomes increasingly attractive as the temperature decreases, eventually violating thermodynamic 
stability criteria. However, as polymer concentration is increased at fixed temperature, the effective 
interaction reverts to mostly repulsive behavior. These issues help illustrate some fundamental 
difficulties encountered when coarse-graining complex systems via effective pair potentials. 



I. INTRODUCTION 



While the computer simulation of single, isolated poly- 
mer chains, either on or off lattice, using a variety of 
conformation sampling algorithms, is nowadays relatively 
routine, for up to L ~ 10 6 monomers or segments 
it is computationally much more demanding to simulate 
polymer solutions or melts, involving large numbers of 
interacting polymer chains. Indeed, if N is the number 
of such chains and L the number of segments per poly- 
mer, then the total number of interacting particles, NL, 
can become very large, particularly so if L is sufficiently 
large for the scaling regime to be reached Q. Under 
these conditions it is tempting to seek a coarse-graining 
procedure to reduce the full segment-level description to 
a model involving only the center of mass (CM) or the 
central monomer of each chain, thus reducing the initial 
AL-body problem to a A-body problem. This is formally 
achieved by tracing out the individual monomer degrees 
of freedom, i.e. by averaging over polymer conforma- 
tions for fixed positions of the CM's or central monomers 
of interacting polymer coils, taking into account the ap- 
propriate Boltzmann weights. This idea goes back to 
Flory and Krigbaum [3j. They predicted that the ef- 
fective repulsive interaction between the CM's of linear 
self-avoiding walk (SAW) polymers should diverge with 
molecular weight at full overlap, i.e. when the CM's of 
the two coils coincide. It was first realized by Grosberg 
et al. that in fact the pair potential between CM's 
remains finite in the scaling limit L —> oo, and of the or- 
der of a few fcgT, reflecting the purely entropic origin of 
the effective interaction. The effective CM pair potential 
between two isolated SAW polymer coils was explicitly 
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evaluated by Monte Carlo (MC) simulations of on and 
off-lattice models , and by renormalization group 
(RG) calculations These studies show that the ef- 
fective pair potential is purely repulsive, with an overlap 
(zero separation) value of about 2k B T and a range of the 
order of the radius of gyration R g . The zero concentra- 
tion pair potential is reasonably well represented by a 
single Gaussian of width R g . 

More recently the effective pair potential was deter- 
mined at finite polymer concentration by a combina- 
tion of monomer-level MC simulations of lattice SAW 
polymers, and an inversion technique based on integral 
equations for pair distribution functions in simple liq- 
uids 0- The resulting pair potentials depend now mod- 
erately on polymer concentration 0, llCf . but they re- 
main essentially repulsive and of range R g . They have 
been put to good use to reproduce the interfacial tension 
of semi-dilute polymer solutions near hard walls or col- 
loidal spheres [ 1 1| . and to determine the depletion inter- 
action between colloidal particles induced by interacting 
(rather than ideal) polymer coils ^2 an d the resulting, 
depletion-induced phase diagram of colloidal dispersions 

m 

The present paper reports an extension of the above 
inversion strategy to the case of dilute and semi-dilute so- 
lutions of interacting linear polymers in solvents of vari- 
able quality, spanning the range between good solvent 
conditions, modeled by the SAW, and poor solvent con- 
ditions, where the coils contract to avoid contact with 
the solvent. This generalization is achieved by adding 
a finite attractive interaction between nearest-neighbor 
monomers of the same or different chains, while main- 
taining infinite repulsion between overlapping monomers. 
The finite nearest-neighbor attraction introduces an en- 
ergy scale, and hence a temperature dependence of the 
effective pair potential. Good solvent conditions cor- 
respond to the infinite temperature limit, leading back 
to the SAW model considered earlier, while increasingly 
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poor solvent conditions are mimicked by enhancing the 
attraction between monomers, or equivalently, decreasing 
the temperature. 

The dependence of the properties of a single polymer 
coil, like its radius of gyration R g , on temperature, and 
in particular around and below the 9 temperature where 
the coil to globule transition takes place, has been stud- 
ied by extensive MC simulations of very long chains (L 
up to 10 6 ) by Grassberger and collaborators [Tj|- Results 
for the temperature dependence of the effective pair po- 
tential for an off-lattice model of two polymer coils (i.e. 
in the infinite dilution limit) were reported by Dauten- 
hahn and Hall 6]. These authors met with increasing 
MC sampling difficulties as the temperature is decreased 
towards the 9 point. The present paper improves on their 
results by considering a lattice model, and using an over- 
lapping distribution method 0, 0] to sample rare con- 
formations which become increasingly important as the 
temperature decreases. The main objective of this paper 
is to consider, for the first time, the effective pair poten- 
tial between the CM's of polymer coils in poor solvent at 
finite polymer concentration, i.e. to investigate both the 
temperature and concentration dependence of the CM 
pair potential, in particular in the vicinity of the 9 point. 
The polymer model studied in this work is defined in Sec- 
tion [n] The zero concentration limit is first considered 
in Section lTm while finite polymer concentrations are ex- 
amined next in Section llVI Some considerations on the 
use of effective interactions in relation to thermodynamic 
stability are presented in Section and conclusions are 
drawn in Section lVll 



II. MODEL AND BASIC PHENOMENOLOGY 

For all the calculations in the present work, a lat- 
tice model defined on the simple cubic lattice (with a 
coordination number of six) was used |17|. A polymer 
chain of length L is represented as a connected sequence 
of L lattice sites. The monomers, defined as occupied 
sites, interact with each other via excluded volume repul- 
sion, preventing two of them from occupying the same 
lattice position, and lattice nearest-neighbor attraction 
e < between non-bonded pairs of monomers. For con- 
venience, the temperature T will be expressed in units of 
—e/kB, or equivalently we will set e = — f and ks = 1. 
With this convention, the usual (3 appearing in the Boltz- 
mann factor is simply the inverse of the temperature, 
i.e. (3 = l/T. In the following, we shall call "contact" 
any nearest-neighbor pair of non-bonded monomers. We 
use "intramolecular" if both monomers belong the same 
chain, and "intermolecular" otherwise. When a system 
consisting of N chains of L segments is considered on a 
portion of cubic lattice of M sites, the monomer packing 
fraction is equal to the fraction of lattice sites occupied by 
polymer segments, c = NL/M, while the concentration 
of polymer chains is p = c/L = N/M. 

This simple model has been extensively studied in 
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FIG. 1: Schematic temperature-inverse length diagram show- 
ing the three regimes of polymer phenomenology and the be- 
havior of the temperatures T$, Tb, Jstabj and T c (defined in 
the text) related to the characterization of the 8 regime. 



order to investigate the properties of isolated polymer 
chains and of polymer solutions. In the phenomenology 
of these systems, three different domains, sketched on 
Fig-ffl are usually defined, corresponding to different be- 
haviors of an isolated chain. 

In the good solvent regime, corresponding in the 
present model to high temperatures, the properties of the 
system are essentially determined by the entropic effects 
originating in connectivity and excluded volume interac- 
tions. An isolated chain takes swollen coil configurations 
with a radius of gyration R g known to scale as L v , where 
v ~ 0.588 is the Flory exponent in three dimensions [Isj . 
This regime has been extensively studied in the case of 
the SAW to which the pre sent model reduces in the infi- 
nite temperature limit [l9| . 

In the bad or poor solvent regime, corresponding here 
to low temperatures, entropic effects are overwhelmed 
by the attractive monomer-monomer interactions and a 
chain molecule collapses into a compact globule with a 
radius of gyration scaling like L 1 ^ 3 . 

Between these two extremes lies the 9 regime, where 
energetic and entropic effects compensate each other. 
In this regime, in the limit of infinite length, a chain 
has the scaling properties of an ideal chain, for instance 
R g cx L 1 / 2 . The characterization of this domain has 
been the focus of many studies using the present lattice 
model [13, 0, H2, , with in particular extensive work 
by Grassberger and coworkers [lj, |2J, |25| . 

Of key importance is the so-called 9 temperature, Tg, 
at which the tricritical coil-to-globule transition takes 
place for an isolated chain in the L — * oo limit. It corre- 
sponds to a real thermodynamic singularity and its best 
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estimate is presently [24| 

T e = 3.717 ± 0.003, i.e. [3 e = 0.2690 ± 0.0002. 

Obviously, simulations can only deal with chains of finite 
length, for which the tricritical singularity is rounded off 
by significant finite size corrections. Tg is thus not di- 
rectly accessible and it is common practice, in simulations 
as well as in experiments, to obtain it by extrapolation 
of related, length-dependent characteristic temperatures. 

The most widely used of such temperatures is Boyle's 
temperature, T B (L), defined as the temperature at which 
the osmotic second virial coefficient of chains of length L 
vanishes. From its knowledge, one can compute the 9 
temperature from the limit 

lim T B (L)=T e . (1) 

L — >OG 

This route has been followed in particular in Ref. [Tij] . 
where it was found that Tb{L) is always greater than Tg 
and decreases with increasing L. 

An alternative possibility is to use the demixing crit- 
ical temperature of the polymer solution. This marks 
the limit between a high temperature regime where the 
chains in solution form a homogeneous fluid, and a low 
temperature domain where the chains start to coagulate 
and the fluid demixes into polymer-rich and polymer- 
poor phases. In the lattice model used here, where the 
solvent is taken into account implicitly through the in- 
troduction of constant attractive monomer-monomer in- 
teractions and variations of the temperature, it is simply 
the critical temperature T C {L) of the liquid-gas transition 
of the system. As for Tb(L), one has 

lim T C {L) = Tg, (2) 

L — >oo 

and it has been found in Refs. |2^ and j2f| that T C (L) is 
always smaller than Tg and increases with increasing L. 

We anticipate on the following sections and introduce 
here a third characteristic temperature called the "sta- 
bility" temperature and denoted T sta h{L). It is associ- 
ated with a fundamental breakdown in the statistical- 
mechanical treatment of the coarse-grained effective pair 
interaction U2<V) between the CM's of two polymer chains 
at low temperature. Singularities of a similar nature have 
been encountered in various models for soft matter sys- 
tems. For instance, Baxter's sticky sphere model |26] dis- 
plays anomalous clustering towards a closed-packed crys- 
talline phase as a result of breakdown of thermodynamic 
stability [2711. In the case of polymers, the Domb- Joyce 
model |28| shows self-trapping behavior, i.e. a chain re- 
mains of finite extent in the infinite length limit, when a 
negative energy is attributed to self-crossings 29]. 

The relation between T sta b(i) and the 9 point is un- 
clear. For any length L, T sta ,b(L) will be shown to be 
smaller than T B (L), the equality being only achieved 
when the pair interaction is zero. For the specific length 
L = 100 studied in detail in this paper, it is found to 
fall between Tb(L) and Tg. Because it corresponds to 



a singularity in the statistical-mechanical treatment of 
the system, it is tempting to assume, as it was done in 
Ref. |33, that 

lim T stab (£) = Tg, (3) 

L — >oo 

and Fig. ^ has been drawn under this assumption. In- 
terestingly, this would imply that, in the infinite length 
limit, since Tg(L — > oo) = T sta b(i - * oo) = Tg, the effec- 
tive interaction between two isolated polymer chains at 
the 9 point vanishes identically for all distances, a behav- 
ior which is trivially found in the case of the Domb- Joyce 
model 28], when a positive constant energy penalty is 
counted for each polymer crossing (note that in this 
model, the temperature scale is reversed since the SAW 
is obtained at zero temperature and ideal chain behavior 
is found at infinite temperature). 

III. ZERO DENSITY LIMIT 
A. Simulation methodology 

In the zero density limit, the effective interaction po- 
tential vi (r) between two polymer chains is equal to the 
difference between the free energy of the two chains with 
their centers of mass constrained to stay at a fixed dis- 
tance r, and the free energy of the same chains infinitely 
far apart. It can be expressed as follows. 

Consider two polymer chains labeled A and B. When 
these chains have conformations T A and T b respectively, 
with the vector tab joining their centers of mass, the 
energy of the pair is 

H(r AB ;T A ,T B ) = iJ intra (r A ) 

+ #intra(r B ) + H ixAer (r AB ; T A , V B ) , (4) 

where the intra- and intermolecular parts (including the 
hard-core interactions) have been separated. Introducing 
the intermolecular Boltzmann weight 

W(\r AB \;T A ,r B ) = exp[-/3H inteI (r AB ;T A ,r B )], (5) 

the effective pair potential is given by |3l| 

pv 2 (r) = -ln(W(\v AB \ =r;T A ,r B )), (6) 

where the brackets denote an average over the probability 
distribution of two isolated chains, which is the square of 
the probability distribution V of a single chain, i.e. 

{W{\v AB \=r;T A ,T B )) = 

V(T A )V(T B )W(\r AB \=r;T A ,T B ) (7) 

TaXb 

with 

V{T) = cxp [-/JiTintraCT)] j £ exp [-/? H intTai (T)} . 

(8) 
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This result provides us with a direct means to com- 
pute the effective pair potential between two chains by 
Monte Carlo simulations. We sample configurations of 
two independent chains using the pivot algorithm |32l | 
and standard Metropolis acceptance rules. The latter 
ensure that the chain conformations are generated ac- 
cording to the probability distribution V . After every 
1000 pivot moves for each chain, we calculate the inter- 
molecular Boltzmann weight (JSJ as a function of the CM 
distance, by moving the polymers towards each other, 
while checking for overlap and counting intermolecular 
contacts. Eventually, fivi(r) is obtained by performing 
the unweighted average of VFdr^s = t\Ya^b) on the 
sample considered: 

JV(V) 

J2W(\r AB \=r;r A ,r B ) 
«r) = -ln^ W) , (9) 

where Af(r) is the number of two chain configurations 
with CM's at distance r sampled during the simulation. 

Using this algorithm, we have been able to compute 
with good statistical accuracy effective pair potentials for 
various chain lengths (L — 100, 200 and 500), provided 
(3 < 0.2. 

Unfortunately, this simple, direct method turns out to 
be inadequate at lower temperatures, where the results 
for the effective potential at short distances are found 
to fluctuate strongly with the sample considered, with 
no significant reduction of the corresponding statistical 
uncertainties when the size of this sample is increased. 
The origin of this problem, already seen in the analo- 
gous off-lattice calculation of Dautenhahn and Hall 0, 
has been carefully analyzed by Grassberger and Heg- 
ger in their calculation of the second virial coefficient 
of the same lattice model as the one studied here [lij . 
At low temperature, the polymer chains start to collapse 
significantly and are thus rather compact. This means 
that at short distances, most two-chain configurations 
obtained with the present direct algorithm will overlap 
and the corresponding intermolecular Boltzmann weights 
are then identically zero. But, in the rare cases where the 
two chains do not overlap, many contacts will usually be 
formed, leading to very low negative intermolecular en- 
ergies and thus huge intermolecular Boltzmann weights 
contributing to the pair potential. The need to average 
such unevenly distributed numbers of very different am- 
plitudes gives rise to the large observed fluctuations and 
renders the direct approach useless. 

Hence more elaborate algorithms must be used for low 
temperatures. In order to choose these new tools, it is 
helpful to recognize the strong similarity between the pre- 
vious scheme and Widom's particle insertion method for 
the computation of the chemical potential |33|. Both 
methods indeed involve the averaging of a Boltzmann 
weight over the equilibrium distribution of some unper- 
turbed system, and in fact the observed breakdown of the 
direct calculation of the effective pair potential at low 



temperature parallels that of Widom's method at high 
fluid densities. In the latter case, an efficient approach 
to solve the problem has been devised by Shing and Cub- 
bins [34|], which belongs to the general class of overlap- 
ping distribution methods first introduced by Bennett to 
compute free-energy differences flH IT^ | . We have thus 
implemented such a method which is described in detail 
in the Appendix. 

Using this histogram method, we have been able to 
extend our calculation of the effective pair potentials be- 
tween chains of length L = 100 to lower temperatures, i.e. 
up to [3 — 0.3. However, the method breaks down when 
applied to longer chains, for mainly two reasons. Firstly, 
the acceptance ratio of the elementary Monte Carlo move 
becomes extremely small, leading to serious ergodicity 
problems in the simulations. Secondly, except for small 
separations corresponding to nearly complete overlap of 
the chains, the two computed histograms do not overlap 
at all, rendering the estimation of the required free energy 
difference quite problematic. The origin of this problem 
can easily be understood on a qualitative level. In the 
original system, the chains attract each other strongly 
and lower their intermolecular energy at low temperature 
by elongating towards each other along the axis joining 
their CM's: there is thus a sizeable volume where the two 
chains overlap with the consequent formation of contacts 
in large numbers |35| . On the contrary, in the reference 
system, only unfavorable excluded volume interactions 
exist between the two partly collapsed chains which are 
thus compressed along the axis joining their CM's to limit 
their overlap : contacts are then formed in small num- 
bers only at the surface separating the two segregated 
chains. This significant separation in the number of con- 
tacts does not occur (or at least is less pronounced for 
the lengths studied here) at complete overlap because 
the system has then to keep the spherical symmetry of 
the isolated chains, thus leaving no possibility for the ad- 
verse changes of the chain shapes seen in the cylindrical 
symmetry. 

Although it would be possible to use alternative meth- 
ods (see e.g. [lj, |23|) to sample even longer chains at 
low temperature, we use mainly L = 100 chains in this 
paper, since we do not expect the qualitative features 
we focus on here to depend strongly on chain length (3(| . 
Of course, quantitative features will vary with and the 
detailed scaling behavior can be quite complex j2J, [25| . 



B. Results 

The effective pair interactions /3«2(r) between two 
polymers of length L = 100 are plotted, for various sol- 
vent qualities, in Fig. [21 The product of the potential by 
Airr 2 is shown as well, because the integral of this quan- 
tity plays an important role for the thermodynamics of 
the system [3(J, |33 , and to emphasize the features of the 
potential at large distances. In this figure and in all the 
following figures representing effective pair potentials, all 
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FIG. 2: Zero density effective pair potentials at different inverse temperatures for polymers of length L = 100 on the simple 
cubic lattice. At each temperature, the :r-axis is scaled with the corresponding R g (0), which is shown in Fig. 01 
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FIG. 3: Radii of gyration of isolated polymer chains of differ- 
ent lengths L on the simple cubic lattice as functions of the 
inverse temperature 0. The data are scaled with the infinite 
temperature values, R a (L = 100) = 6.44, R g (L = 200) = 
9.76, and R g {L = 500) = 16.84. 



distances have been scaled at each temperature by the ra- 
dius of gyration of the isolated chain at this temperature, 
which is plotted in Fig.|21for completeness. 

One can distinguish two temperature domains. In 
Ref. , where the infinite temperature limit of the model 
was studied, 0V2(r) was found to be purely repulsive with 
a roughly Gaussian shape centered on r = 0. Here we 
find that it remains so, provided the temperature is high 
enough, i.e. < 0.2. The main effect of lowering the 



temperature is a reduction of the overall amplitude of 
the potential and a slight decrease in its range in terms 
of the normalized distance. 

At lower temperatures, > 0.2, 0V2{r) begins to ex- 
hibit qualitative changes: while the potential retains a 
repulsive Gaussian-like component at short distance, a 
negative, attractive tail appears at large distance. As 
the temperature decreases further, the amplitude of the 
repulsive core decreases, following the trend of the previ- 
ous high temperature domain, whereas the attractive tail 
becomes more and more important, ultimately dominat- 
ing the whole picture, as seen for = 0.3, where 0V2(r) is 
everywhere negative and only a modest repulsion shows 
up for r/Rg < 1. 

A second, less prominent feature appears at low tem- 
perature (0 > 0.2) as well. At r = 0, the effective pair po- 
tential displays a small minimum, which becomes deeper 
when the temperature decreases. This means that full 
overlap is locally stable and that one has to overcome 
a (modest) free-energy barrier to separate two chains in 
this configuration. For chains of length L = 100, it is 
difficult to know if this feature is a generic property of 
the effective pair potential or just a lattice artifact due 
to the shortness of the chains (the width of the minimum 
is indeed of the order of the lattice spacing). However, 
despite our limited ability to investigate longer chains at 
low temperature, we have seen evidence of a similar min- 
imum, but of larger width, for chains of length L = 500 
at = 0.225, suggesting that this is a genuine physical 
effect. 

One can compute various interesting scalar quantities 
from the knowledge of the effective pair potential. Ac- 
cording to Sec. [Hj the second osmotic virial coefficient, 
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FIG. 4: Reduced second osmotic virial coefficient (B2/RI) 
and reduced stability integral at zero density {I2/ R A g ) as func- 
tions of the inverse temperature j3 for polymers of length 
L — 100 on the simple cubic lattice. 



FIG. 5: Zero density effective pair potentials at different in- 
verse temperatures f3 (from top to bottom, for a given length, 
/3 = 0, 0.1, 0.2, and 0.3) for polymers of various lengths L on 
the simple cubic lattice. At each temperature, the x-axis is 
scaled with the corresponding R g {/3). 



given by 

1 f°° 

B*(P) = 7: [1 - exp(-/3 W2 (r))] 4ir 2 rfr, (10) 
z Jo 

is of particular importance. In addition, we introduce the 
"stability" integral, 

1 f°° 

h{!3) = t;J o (3v 2 (r)4wr 2 dr, (11) 

which results from the linearization of the exponential 
in Eq. I|10|l . For systems interacting through density- 
independent potentials, the sign of I2 gives a neces- 
sary condition for the existence of the thermodynamic 
limit and the stability of the system against coalescence 
[ssl l8?il ]. as will be further developed in Sec. E] Ac- 
cordingly, as discussed by anticipation in Sec. ^ a "sta- 
bility" temperature T sta b(i) can be defined, at which 
hiP) vanishes for chains of length L, and below which 
the necessary condition for stability against coalescence 
is violated. The ordering T stah (L) < T B {L) discussed 
above is then an immediate consequence of the fact that 
x > 1 — exp(— x). 

Both quantities, normalized with the cube of the ra- 
dius of gyration to obtain dimensionless quantities, are 
plotted in Fig.^J As was implicitly assumed above, they 
are decreasing functions of (3 and one finds, for the chains 
of length L = 100 studied here, that T B > T sta b ^ 3.8. 

Finally, we address the question of the dependence 
of the zero density effective pair potential on the chain 
length. In Ref. 0, the infinite temperature case has been 
studied in detail over a wide range of polymer lengths 
(L = 100 to 8000). Such an extensive study was not at- 
tempted here and we only considered two rather short 



lengths (L = 200 and 500) in addition to the previous 
data for L = 100. All these data, plotted in Fig. [5] with 
the one previously described, have been computed with 
the direct method, leading to strong statistical uncertain- 
ties at the lowest temperature corresponding to (3 = 0.3. 

To interpret these potentials one should keep in mind 
that R g scales differently with temperature for the dif- 
ferent lengths, as illustrated in Fig. |3] Also, Tb and T c 
are different for each length. Following Refs. 0] and 
PH. in which these temperatures have been computed 
with great care, for L = 100: T c ~ 3.10; T B ~ 3.86; for 
L = 200: T c ~ 3.33; T B ^ 3.81; for L = 500: T c ~ 3.41; 
T B ~ 3.77. 

A few qualitative observations can be made, at least 
for (3 < 0.2, where the statistics are good enough. With 
the normalization chosen for the x-axis, no intersection 
between potentials corresponding to different values of L 
at the same temperature is found. In Ref. this fea- 
ture was already found for (3 = 0, combined with the fact 
that a larger L leads to a less repulsive potential. We 
recover this result here and the same qualitative behav- 
ior is found for (3 — 0.1, but the potentials for different 
values of L are closer. This behavior changes for (3 = 0.2, 
where we now see that the larger L, the more repulsive 
the potential. Globally, we thus find that in the high- 
temperature domain, the decrease in the amplitude of 
the repulsive effective potential is slower for large values 
of L. 

It is difficult to draw conclusions from the data at (3 — 
0.3, but it looks like the larger L, the deeper the attrac- 
tive potential. Indeed, for large CM distances, at which 
the statistical inaccuracies in the simulation results are 
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expected to be modest, the potential seems to be more 
negative when L is large; the same holds at complete 
overlap: trying to estimate ftv^{Q) for L = 500 with the 
overlapping distribution method, we find ftx>2{0) ~ —1.1, 
to be compared to /3i>2(0) ~ —0.15 for L = 100. 

All these results, if applied to Eq. (EH , are consis- 
tent with the findings of Grassberger and Hegger for the 
variation of the reduced second osmotic virial coefficient, 
B 2 /Rg, with temperature for various lengths (see Fig. 16 
in Ref. ^3). Indeed, they found that the larger L, the 
flatter this quantity is in the high temperature regime, 
the more abrupt the downward bend of the curve when 
approaching Tg(L) and the faster the divergence towards 
— oo in the poor solvent regime. 

IV. FINITE DENSITIES 
A. Methodology 

Having derived the effective potential between two iso- 
lated polymers, we now turn to polymers at finite density. 
For this, we follow the route proposed in previous work 
0, in which an effective pair potential was constructed 
to exactly reproduce the two-body CM correlations of 
the full underlying many-body system. In fact, it can be 
proven for a wide variety of systems that for any given 
pair distribution function g(r) at given inverse tempera- 
ture ft and density p, there exists a corresponding unique 
two-body pair potential v(r) which reproduces g(r) irre- 
spective of the underlying many-body interactions in the 
system |40j . Of course, g(r) depends on density and tem- 
perature and contains contributions not only from the 
bare pair potential U2(r), but also from the three- and 
more-body terms. As a consequence, the effective pair 
interaction will also be state dependent and a new effec- 
tive potential, hereafter denoted v(r;p,ft), must be cal- 
culated for each density and temperature. However this 
inversion approach says nothing about a possible volume 
term v\(p,ft), in the coarse-grained total potential en- 
ergy, which contributes to the equation of state, but not 
directly to the pair-correlations 0, • Of course the 
volume terms may still contribute indirectly, for example 
when they induce phase transitions. 

The inversion procedure, using g(r) to extract v(r), 
is well known and has been studied extensively in the 
field of simple fluids |43l l44j. We invert g(r) using the 
hypernetted chain (HNC) closure, 

g(r) = exp(-ftv(r) + g(r) - c{r) - 1), (12) 

of the Ornstein-Zernike equation While the simple 
HNC inversion procedure would be inadequate for dense 
fluids of hard core particles, where more sophisticated 
closures or iterative procedures are required El > we 
are able to demonstrate the consistency of the HNC in- 
version in the present case. 

To compute the necessary structural information we 
have performed canonical Monte Carlo simulations of 



polymer solutions. We have studied chains of length 
L = 100 in a cubic box of size M = 100 3 with peri- 
odic boundary conditions, varying the number of poly- 
mers from N = 400 to N = 3200 (c = 0.04 to 0.32). 
Four temperatures have been considered corresponding 
to ft = 0, 0.1, 0.2, and 0.3. Note that the lowest tem- 
perature, T = 3.33, is slightly larger than the critical 
temperature of the system, T C (L = 100) ~ 3.1 according 
to Ref. |25|. thereby avoiding concerns with the possible 
two-phase behavior of the system. To sample the config- 
uration space of the system we have used standard tech- 
niques of polymer simulations: pivot moves |32| . trans- 
lation moves, and, for high densities where the previ- 
ous moves become inefficient, configurational bias Monte 
Carlo (CBMC) moves, in which an extremity or part of 
the interior of a chain are regrown |l6j . 

In the course of the simulations, the CM of each poly- 
mer was tracked in order to construct the CM radial pair 
distribution function g(r; ft, p). The latter is only known 
up to a cutoff radius r c , which corresponds to half the size 
of the simulation box. For the inversion, we need g(r) for 
all r, so we employ the following iterative scheme to ex- 
tend g(r). As an initial step, we set g(r) = 1 for r > r c 
and calculate the corresponding v(r) by inversion. We 
then set v(r) = for r > r c and determine the corre- 
sponding g(r) for < r < oo by a regular HNC calcu- 
lation, using a simple iterative procedure. The g{r) for 
r < r c is then replaced by the measured g(r), and the new 
v(r) is calculated. This is again set to zero for r > r c , 
and the process is repeated until convergence. In fact, 
because of the finite box-size, the inversion process is un- 
derdetermined, and our ansatz that v(r) = for r > r c is 
needed to find a unique solution 10] . This is not unrea- 
sonable since we do not expect the interactions between 
the polymer coils to be significant beyond distances a few 
times the radius of gyration. However, to make sure that 
this is actually the case, we found that relatively large 
simulation boxes were needed, with a lattice size of up to 
10-15i? 9 . This is particularly important at high density, 
where the inverted potential becomes longer ranged and 
more sensitive to small changes in the radial distribution 
function g(r). In all our inversions, we checked explicitly 
that v(r) becomes effectively zero before r = r c , confirm- 
ing our initial ansatz. 



B. Results 

Before presenting the results for the effective pair po- 
tentials, we first discuss briefly the behavior of the ra- 
dius of gyration with density and temperature, which is 
an important measure of the physical properties of the 
polymer. The corresponding simulation data are plotted 
in Fig. [fjj As found in similar previous work [45|. for 
a given temperature, R g decreases when the density of 
polymer increases if this temperature falls into the good 
solvent regime, whereas R g increases if this temperature 
is located below the 9 regime (Fig. |SJi) . When R g is plot- 
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FIG. 6: Density- and temperature-dependence of the radius of gyration of polymers of length L = 100 on the simple cubic 
lattice, (a) Radii of gyration as functions of the monomer density c for various inverse temperatures (3. (b) Radii of gyration 
as functions of the inverse temperature (3 for various monomer densities c. 



ted as a function of temperature for various densities, 
this results in the existence of a regime, located in the 
9 region, where R g is nearly independent of c, as can be 
seen in FiglSJs, where all the curves are seen to converge 
around [3 ~ 0.25. 

The effective temperature- and density-dependent pair 
potentials are plotted in Figs. and |S] for two representa- 
tive temperatures, corresponding to j3 = 0.1 and (3 = 0.3 
respectively. From these data, it is clear that the two 
temperature domains identified in the zero density study 
of Sec. IIIII display obvious distinctive features in terms of 
the density dependence of the effective pair potentials at 
fixed temperature. 

When j3 = 0.1 (Fig. 0), i.e. in the high temperature 
regime, a moderate density dependence is found. Start- 
ing from the purely repulsive, nearly Gaussian shape of 
the potential at zero density and increasing the density, 
one finds a slight increase in the amplitude and in the 
range of the repulsive potential for small separations of 
the polymer CM's and the progressive appearance of a 
weak negative attractive tail at large distance. This be- 
havior is similar to the one found in the infinite tem- 
perature limit studied in Ref. 0, with the only minor 
difference that we do not see a decrease of the ampli- 
tude of the potential at short distances for the highest 
concentrations . 

The situation is completely different when [3 = 0.3 
(Fig. |SJ. Here the density dependence of the effective 
pair potential is very important and leads to significant 
qualitative changes in the shape of the potential: whereas 
at low density, v(r) is negative and essentially attractive, 
repulsion between centers of mass progressively builds 
up when the density increases, so that an almost purely 
repulsive, Gaussian shaped v(r) is eventually observed 



at the largest density. For the good solvent case a clear 
link has been found between the density dependence and 
the strength of the 3- and higher body interactions [4(1 
|47| . Presuming that the same link can be made here, 
this would imply that the many-body interactions are 
relatively more important at lower temperatures. 

The virial coefficient has no clear interpretation for 
potentials obtained at finite global densities, whereas the 
"stability" integral does still define a lower limit to the 
existence of a thermodynamic limit of the coarse-grained 
system. As can be directly seen from Fig. |SJd, at all den- 
sities, except the largest (c = 0.32), the stability integral 
of the computed potentials is negative, leading to serious 
consistency issues to be discussed in the next section. 



V. CONSIDERATIONS OF THERMODYNAMIC 
STABILITY 



The potentials computed in this work raise a num- 
ber of interesting conceptual issues with regard to the 
use of effective potentials in coarse-grained descriptions 
of material properties. The first question concerns ther- 
modynamic stability |38(: do the calculated pair poten- 
tials generate a valid thermodynamic limit ? Secondly, 
there is a question of representability |47| : for a given 
state point, how well can the properties of the underly- 
ing polymer system be coherently represented by a single 
coarse-grained effective pair potential ? 
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A. Thermodynamic stability of effective potentials 

We briefly repeat the criteria for the existence of a 
thermodynamic limit described in Ruelle's classic book 
|38j and valid for state independent interactions. Con- 
sider a system of N particles in a volume V. If the total 
interaction energy Vn, which can be built from pair and 
higher order terms, satisfies, for all N > and for all 
configurations {r^} in the configuration space R N , the 
inequality 

V N (r u ...,r N )>-BN (13) 



with B > 0, then, according to definition 3.2.1. in 
Ref. [3^, the system is stable: the grand partition func- 
tion converges, and there is a well-defined thermody- 
namic limit. Potentials that do not satisfy this criterion 
are termed catastrophic by Ruelle. Specializing to pair 
potentials, i.e. 

l4 2) (r 1 ,...,r A r)= Yl (14) 

l<i<j<N 

condition (|13|l leads to the following necessary (but not 
sufficient) condition for the existence of a thermodynamic 
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limit: 

v(0) = J v 2 (r)4:iTr 2 dr > 0, (15) 

which, following Eq. Ijllfl . can be rewritten I 2 > 0, hence 
the name "stability" integral for 12- If Eq. (|15fl is not 
satisfied, then condition i|13|l can be violated for con- 
figurations of a homogeneous fluid. As N — > 00, the 
free energy grows super-extensively - the system has no 
thermodynamic limit - and particles coalesce to form a 
dense cluster. Condition i|15|) is necessary, but not suffi- 
cient, for the existence of a thermodynamic limit, since 
one can also construct potentials for which -0(0) > 0, but 
where the system in a microscopically inhomogeneous 
(typically, crystalline) state is unstable to coalescence. 
See Refs. [38| and [3(| for explicit examples. 

Effective potentials with hard cores, such as those used 
to describe simple atomic and molecular materials, can 
be easily shown to satisfy the criterion l|13l) . In contrast, 
the potentials describing the effective interactions be- 
tween the CM's of polymers studied in the present work 
do not have a hard core, leading to the possibility that 
two or more "effective" particles occupy the same posi- 
tion in space. The existence of a thermodynamic limit, 
where the free energy per particle is bounded, is therefore 
much more subtle. 

Although the potentials calculated in this work are 
state dependent and technically only relevant at the den- 
sity for which they have been obtained, it is nevertheless 
interesting to consider what would happen should one of 
them be used, independently of the density, to describe a 
system of N > 2 particles in a volume V. To this end, we 
first consider the zero density potentials shown in Fig. [21 

When the underlying SAW polymer system has no 
nearest-neighbor monomer attractions, the effective pair 
potentials were shown to be positive 0, H3 , so that the 
criterion (|T3|l is obeyed. However, as Fig.|21demonstrates, 
the introduction of nearest-neighbor attractions leads to 
effective potentials which are no longer positive definite. 
As the temperature is lowered, the potential grows more 
and more attractive, until finally it violates the Ruelle cri- 
terion and becomes catastrophic. This can be diagnosed 
in Fig. 01 where 1% becomes negative at low temperature 
and the necessary condition (|15|) is thus violated. In prin- 
ciple, the stability limit of the system should be traced by 
looking for a state point at which the potential leads to 
coalescence into microscopically inhomogeneous states, a 
rather difficult task in general. Here we will use a sim- 
pler, approximate criterion, namely Eq. i|15|) . which is a 
necessary, but not sufficient condition for stability. We 
will make the heuristic assumption that the true stability 
limit is not far removed from the simpler one we employ 
[48| . For a given potential then, Eq. 1|15|) allows us to 
define the "stability" temperature T sta b below which the 
effective pair potential becomes catastrophic. This ex- 
plains our nomenclature in Sees. ITU and IIIII 

The potentials derived at zero density violate Eq. i|15|) 
below a temperature T sta b ~ 3.8; the potentials are ex- 



pected to become unstable to inhomogeneous coalescence 
at some temperature above that. If one were to use po- 
tentials derived for T < T sta b at a finite density, the 
system would be catastrophic. 

This naturally leads to the next question: what about 
the potentials we derived at finite density ? There, 
as shown in Fig. Q3J the pair potential can also violate 
Eq. (|15ll . even though the underlying polymer system is 
stable, both to coalescence and phase separation, since 
T > T c for our L = 100 SAW lattice chains. For large 
enough densities c, the potentials again satisfy Eq. (|15|l . 

Serious difficulties are thus emerging with the present 
coarse-graining procedure: in certain portions of the 
temperature-density plane, we replace a complex, but 
perfectly well-behaved, polymer system by a simple soft 
colloid fluid with pathological thermodynamics! This 
also raises a certain number of formal concerns. The 
first one is about the unicity of the mapping between 
g(r) and v(r) invoked in Sec. IIVI and demonstrated in 
Ref. ^(j- The theorems rely on well-defined statistical 
ensembles and their derivations are no longer valid when 
a catastrophic potential is the outcome of their applica- 
tion. Similarly, the lack of a well defined ensemble also 
leads to questions about the validity of the liquid state 
theory, including the Ornstein-Zernike equation, etc. 

However, in spite of these formal difficulties, one could 
take a purely pragmatic attitude: Since it is numerically 
possible to extract effective pair potentials using the pro- 
posed procedure, why not just ignore all the previous 
concerns and see if these potentials can be of any practi- 
cal use ? To investigate this point, MC simulations of a 
soft colloid system interacting through the catastrophic 
effective pair potential obtained for polymers at /3 = 0.3 
and c = 0.04 were performed, revealing a very interesting 
behavior. 

We indeed find that the behavior of the soft colloid 
fluid depends strongly on the configuration from which 
the simulation is initiated. When an initial configuration 
is set up with all soft particles gathered in a single dense 
cluster, the system coalesces: the particles stay together 
indefinitely (after a few million MC steps per particle, 
not a single one has been seen to escape the cluster) , and 
the potential energy per particle e appears unbounded, 
increasing in absolute value with increasing N; more pre- 
cisely, we find e ~ -24.6 for N = 400 and e ~ -202 for 
N = 3200, which would give e(N) ~ -0.062A. This is 
exactly the type of behavior expected from a system in- 
teracting through a catastrophic potential. The structure 
of the clusters also shows distinctive features, as can be 
seen in Fig. EI where pair distribution functions f(r) are 
plotted for the previous cluster sizes. For N = 3200, /(r) 
displays three peaks, located at r/R g = 0, r/R g ~ 1, and 
r/Rg ~ 1.5. Quite evidently, the first two peaks origi- 
nate in the formation of groups of superimposed particles 
separated by the distance corresponding to the interac- 
tion potential minimum. This is indeed a very efficient 
way for the system to lower its energy, since the energy 
cost for overlapping particles is modest compared to the 
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FIG. 9: Pair distribution functions of two collapsed soft col- 
loid systems, with N = 400 and N = 3200 particles respec- 
tively. In both cases, the space integral of f(r) has been 
normalized to 1 to allow a more significant comparison. The 
interaction potential is the effective pair potential obtained 
for a polymer system at /3 = 0.3 and c = 0.04. 



stabilization of pairs of particles at distance r ~ R g . As 
for the third peak, it is very likely it has to be associ- 
ated with the position of the second nearest neighbors. 
In comparison, for TV = 400, f(r) is rather featureless 
and longer-ranged. This change in shape of f(r) is again 
easily understood from the clustering mechanism. For 
small N, smaller groups of superimposed particles can 
be formed, creating less deep energy wells at distance 
r ~ R g ; the resulting clusters are thus more diffuse and 
less structured than for large N. 

But, if the initial conditions correspond to a homo- 
geneous distribution of particles in the simulation box, 
the fluid is found to remain homogeneous over the entire 
time of our simulations (up to 45. 10 6 MC steps per par- 
ticle for 400 particles), suggesting that the homogeneous 
fluid is metastable with respect to catastrophic coales- 
cence. In this case, the potential energy per particle e is, 
as expected, independent of the size of the system, i.e. 
e ~ —0.176 for systems containing 400 or 3200 particles. 

Of course, the original polymer system, which is per- 
fectly stable, does not show such a dependence of its 
behavior on the initial conditions. This can be easily 
checked by preparing an initial configuration consisting 
of a dense cluster of polymer chains: the cluster "evap- 
orates" very rapidly, and the density becomes uniform 
again within a few thousand MC steps per chain only. 

In Fig.[23 we compare the g(r) of the metastable fluid 
phase of the catastrophic soft colloid system to the g{r) 
of the CM's of the underlying polymer system at c = 0.04 
and f3 = 0.3. The agreement is excellent! This suggests 
that the HNC closure used in the inversion procedure, 
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r/Rg 

FIG. 10: Comparison of the pair distribution functions ob- 
tained from simulations of a polymer system (symbols) and of 
the corresponding soft colloid fluid in its homogeneous phase 
(solid line). Here, (3 = 0.3, N = 400, L = 100, M = 100 3 , 
and hence c = 0.04. 

a closure which has been shown to be very accurate for 
bounded stable potentials still works, in spite 

of the fact that the potential it produces is catastrophic 
p^|. The situation is akin to that of the Percus-Yevick 
solution for the structure of Baxter's sticky sphere model 
pfij , which is useful for describing hard colloidal systems 
with very short ranged attractive interactions in spite 
of the fact that the underlying model system is actually 
catastrophic j^]. In both cases, the problem is most 
likely circumvented by the approximate character of the 
chosen closure. 



B. Representability problems for effective 
potentials 

The considerations of the previous section, where the 
coarse-graining procedure led to catastrophic potentials, 
are examples of the more general problems of transfer- 
ability and representability of effective potentials used to 
describe complex systems, issues discussed in more detail 
in a recent review |47| . 

Transferability problems occur when an effective po- 
tential derived at one state point is not applicable at a 
different state point. In essence all the derived effective 
potentials in this paper suffer from transferability prob- 
lems, since they vary with density and with temperature. 
However, for high temperatures, the density dependence 
is not that strong, so that the transferability problems 
are not as important as they appear to be at lower tem- 
peratures. For example, the effective potentials in Fig. [S] 
vary much more rapidly with density, suggesting that one 
must be very careful in using a potential derived at one 



12 



state point as an approximate one for another state point. 

Representability problems occur when effective poten- 
tials derived to reproduce one physical property, do not 
easily describe another physical property |47j . For exam- 
ple, one can calculate the (unique) pair potential v g (r) 
that reproduces the homogeneous structure (described 
by g(r)) of a system interacting via a Hamiltonian with 
pair and triplet interaction terms (similarly to what was 
done in this paper) . If the usual equations for the internal 
energy or the virial pressure, valid for normal pair poten- 
tials, are applied to v g (r), then the latter thermodynamic 
quantities are not correctly reproduced. This was pointed 
out many years ago in the context of simple fluids |50| . 
The same representability problems were already found 
for the athermal polymer case ^3] • They are expected to 
be more pronounced when the density dependence is as 
prominent as in Fig. [5] Moreover, as pointed out in the 
previous section, even more vexing representability prob- 
lems occur, because the derived effective potentials can 
lead to systems without a well-defined thermodynamic 
limit. 

The problems of transferability and representability 
are more important for inhomogeneous systems. For ex- 
ample, at any interface between two phases, it is not clear 
which potential should be used. Similarly, for the tem- 
perature regime that leads to the catastrophic potentials 
of the previous section, the apparent metastability of the 
effective "soft colloid" system allows for the inversions to 
work for a homogeneous system, but this breaks down if 
the simulations are started with certain inhomogeneous 
initial conditions. Using some measure of the local den- 
sity may be a better path to follow. Taking again the 
example of Fig. HO if one were to use a local density de- 
pendence, the effective system could be stable against 
collapse, as the potentials would become more repulsive 
for higher local densities. This would then more closely 
resemble the underlying polymer system. However, pre- 
scriptions for taking local density dependence into ac- 
count which are both accurate and tractable are not yet 
well developed. 

A final issue not yet resolved is the possible role of 
volume terms - contributions to the free-energy which are 
independent of the particular CM configuration. Their 
effects on phase behavior can be subtle (see e.g. |4ll 142^ 1. 
and they may appear in certain coarse-graining schemes 
[5l| . For the case of polymers in a good solvent, they were 
shown to be negligible @, E3 > but that may no longer be 
the case for poorer solvents. 



VI. CONCLUSION 

We have extended previous work on a coarse-grained 
description of polymer solutions in good solvent to in- 
creasingly poor solvent conditions by adding nearest- 
neighbor attractions to the initial lattice SAW model, 
and gradually lowering the temperature. As in the ear- 
lier work |t1 fill |4^|. an effective pair potential v(r) be- 



tween the centers of mass of the linear polymer coils was 
extracted from the Monte Carlo generated CM-CM pair 
distribution function g(r), using the HNC closure, which 
is known to be very accurate in the absence of any hard 
core repulsion. In the infinite dilution limit, the effec- 
tive pair potential V2(r) is close to a Gaussian-shaped 
repulsion, of amplitude /3v2(0) « 2, at high tempera- 
tures (corresponding to the SAW limit), but this ampli- 
tude decreases as the inverse temperature f3 increases; 
an attractive tail develops and at the lowest temperature 
investigated in our MC simulations, the effective pair po- 
tential is entirely attractive, signalling a tendency of the 
system to coalesce. 

At these low temperatures, severe ergodicity problems 
arise in the MC simulations, and, in the zero density 
limit, an overlapping histogram method must be used 
to extract statistically significant results. The problem 
worsens with increasing polymer length L, so that our 
simulations were mostly restricted to L = 100. 

Increasing the polymer concentration at fixed temper- 
ature leads to a "restabilization" of the solution in the 
sense that the effective pair potential exhibits an increas- 
ingly repulsive component as the system is taken from 
the ultra-dilute to the semi-dilute regime. Even at the 
lowest temperature investigated ((3 — 0.3), the effective 
pair potential reverts to an almost exclusively repulsive 
Gaussian-like shape at c = 0.32 (which, for these poly- 
mers, is close to the melt regime). 

The occurrence of strongly attractive and significantly 
state-dependent effective pair potentials between poly- 
mer CM's raises the question of the thermodynamic sta- 
bility of systems of particles interacting via such "catas- 
trophic" potentials, and of the appropriateness of the 
coarse-graining procedure to describe solutions of inter- 
acting polymers, which are intrinsically stable, close to 
the 9 temperature. The observation that the effective 
coarse-grained pair potential v(r) is capable of repro- 
ducing the g(r) (derived from monomer scale simula- 
tions) despite the catastrophic nature of the pair poten- 
tial, points to the possible existence of metastable ho- 
mogeneous states generated by these potentials. These 
metastable states appear to exhibit proper thermody- 
namic extensivity properties, and may be stabilized 
against ultimate coalescence by sufficiently high kinetic 
barriers. If this is indeed the case, the catastrophic ef- 
fective pair potentials may still provide a useful coarse- 
graining tool to describe homogeneous states. Strongly 
inhomogeneous states, generated by coalescence, lead to 
widely different local densities, and hence would require 
the use of effective pair potentials depending on the local 
density. The fact that the effective pair potential tends 
to become more repulsive at higher density might provide 
the stabilizing mechanism for the homogeneous state. 

We plan to extend the present work to examine poly- 
mer CM density profiles near interfaces, and to extract 
the osmotic equation of state of polymer solutions as a 
function of concentration and temperature. 
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APPENDIX A: OVERLAPPING HISTOGRAM 
METHOD FOR THE COMPUTATION OF THE 
EFFECTIVE POTENTIAL BETWEEN TWO 
CHAINS 

We first introduce a reference system consisting of two 
polymer chains similar to the original ones, except that 
the intermolecular interaction does not include an attrac- 
tive nearest-neighbor contribution. The corresponding 
effective pair potentials v^ir) can be efficiently com- 
puted with the direct method, since the intermolecular 
Boltzmann weights can only be either or 1 in this case. 

The overlapping distribution method is then used to 
compute Viir) — v\ (r), the free-energy difference be- 
tween the two original chains, with their CM's con- 
strained to stay at a given distance r, and the two ref- 
erence chains under the same geometrical constraint. To 
do so, we need to know for these systems the distribution 
functions 7r or i g (r, n c ) and 7r rc f(r, n c ), resp ectively, of the 
number n c of intermolecular contacts |52f. which, accord- 
ing to the theory of the overlapping distribution method, 
obey the relation 

(3v 2 (r) ~(3vf(r) = In W r '" c ) + (3en c , (Al) 

7T re f(r, n c ) 

where e = — 1 denotes the nearest- neighbor pair attrac- 



tion. 

These histograms have been computed by sampling in 
canonical Monte Carlo simulations, for both the origi- 
nal and the reference systems, the ensemble of two chain 
configurations with CM's at distance r. The elementary 
move we use is the following. One chain is chosen at 
random and a pivot move is attempted on it. If no self- 
overlap, leading to immediate rejection, occurs for this 
chain, a random position with CM's at distance r is cho- 
sen for the modified chain around the other coil. If in 
this position the two chains overlap, the move is rejected, 
else it is accepted or rejected according to the standard 
Metropolis rules depending on the change in total energy. 

Examples of the accumulated histograms, together 
with their combination through formula IjAlfl , are shown 
in Fig. ^2 Clearly, the explicit dependence of the left- 
hand side of Eq. I|A1J| on n c is found to disappear as 
prescribed by the theory, giving us a good indication of 
the convergence of the method. 

The corresponding zero density pair potentials for 
polymer chains of length L = 100 are shown in Fig. 1121 
where they are compared to the results obtained with the 
direct simulation method. j3 = 0.2 is the largest value at 
which the direct method gives smooth and reproducible 
results and thus has been used as a test case for the 
overlapping histogram approach: as can be seen on the 
figure, the curves obtained with both simulation meth- 
ods are indistinguishable. For larger /3's, the effective 
pair potentials computed with the direct method display 
irregularities at small distances, especially apparent for 
f3 = 0.3, and originating in the sampling difficulties men- 
tioned in Sec. IIIII This is clearly not the case of those ob- 
tained with the histogram approach which are perfectly 
smooth and meet the preceding curves for large enough 
distances at which sampling problems disappear. 
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FIG. 11: Distribution functions for the number n c of intermolecular contacts between two simple cubic lattice chains of length 
L = 100 at /3 = 0.3 and with CM's constrained to stay at distance r — 4.2. 7r or i g (r, n c ) and n le {(r,n c ) refer to the original 
and reference systems respectively, (a) Log plot of 7r rig(r, n c ) and 7r rc f(r, n c ) as functions of n c . (b) Determination of the 
free energy difference V2(r) — V2 el (r) using Eq. 1A1I . As expected from the theory, no explicit dependence of the combination 
In [7Torig(r, n c )/-K val (r, n c )] + (3en c on n c is found. 
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FIG. 12: Comparison of the zero density effective pair poten- 
tials obtained with the overlapping histogram method (sym- 
bols) and with the direct approach (solid lines) at various 
inverse temperatures /3 for polymers of length L = 100 on the 
simple cubic lattice. 



